# -*- coding: utf-8 -*-
"""
Created on Mon Feb  4 07:36:00 2019

@author: Salem
"""

import numpy as np
import numpy.linalg as la


import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


#-------------------------------------
from pathlib import WindowsPath
from pathlib import Path
import os

cwd = Path(os.getcwd())


result_dir = cwd  / "Run Results" 

matrix_dir = result_dir / "HD Results" / "Matrices_WTRC_Affine"

mesh_dir = result_dir  / "MHD Results" / "CH_Affine" / "Meshes"
#-------------------------------------

import Aligning_Meshes as AM


#CH_ICP_distances_Rigid = np.load(result_dir / "CH_ICP_distances_Rigid.csv")
#D_CH_ICPR = 0.5 * (CH_ICP_distances_Rigid + np.transpose(CH_ICP_distances_Rigid))



ICP_order = [17, 18, 19, 23, 24, 25, 20, 21, 22, 14, 15, 16,26,27, 28, 29,0, 1, 2, 3, 4, 5, 6, 7, 8 ,11, 12, 13,
             9, 10, 30, 31, 32]


CH_ICP_distances_Rigid = np.loadtxt(open(result_dir /"CH_ICP_distances_Rigid.csv", "rb"), delimiter=",")
CH_ICP_distances_Rigid[:,:] = CH_ICP_distances_Rigid[ICP_order,:]
CH_ICP_distances_Rigid[:,:] = CH_ICP_distances_Rigid[:,ICP_order]
D_CH_ICPR = 0.5 * (CH_ICP_distances_Rigid + np.transpose(CH_ICP_distances_Rigid))

CH_ICP_distances_Uniform = np.loadtxt(open(result_dir /"CH_ICP_distances_Uniform.csv", "rb"), delimiter=",")
CH_ICP_distances_Uniform[:,:] = CH_ICP_distances_Uniform[ICP_order,:]
CH_ICP_distances_Uniform[:,:] = CH_ICP_distances_Uniform[:,ICP_order]
D_CH_ICPU = 0.5 * (CH_ICP_distances_Uniform + np.transpose(CH_ICP_distances_Uniform))

CH_ICP_distances_Resize = np.loadtxt(open(result_dir /"CH_ICP_distances_Resize.csv", "rb"), delimiter=",")
CH_ICP_distances_Resize[:,:] = CH_ICP_distances_Resize[ICP_order,:]
CH_ICP_distances_Resize[:,:] = CH_ICP_distances_Resize[:,ICP_order]
D_CH_ICPS = 0.5 * (CH_ICP_distances_Resize + np.transpose(CH_ICP_distances_Resize))

CH_ICP_distances_Affine = np.loadtxt(open(result_dir /"CH_ICP_distances_Affine.csv", "rb"), delimiter=",")
CH_ICP_distances_Affine[:,:] = CH_ICP_distances_Affine[ICP_order,:]
CH_ICP_distances_Affine[:,:] = CH_ICP_distances_Affine[:,ICP_order]
D_CH_ICPA = 0.5 * (CH_ICP_distances_Affine + np.transpose(CH_ICP_distances_Affine))

WTRR_ICP_distances_Rigid = np.loadtxt(open(result_dir /"WTRR_ICP_distances_Rigid.csv", "rb"), delimiter=",")
WTRR_ICP_distances_Rigid[:,:] = WTRR_ICP_distances_Rigid[ICP_order,:]
WTRR_ICP_distances_Rigid[:,:] = WTRR_ICP_distances_Rigid[:,ICP_order]
D_WTRR_ICPR = 0.5 * (WTRR_ICP_distances_Rigid + np.transpose(WTRR_ICP_distances_Rigid))

WTRR_ICP_distances_Uniform = np.loadtxt(open(result_dir /"WTRR_ICP_distances_Uniform.csv", "rb"), delimiter=",")
WTRR_ICP_distances_Uniform[:,:] = WTRR_ICP_distances_Uniform[ICP_order,:]
WTRR_ICP_distances_Uniform[:,:] = WTRR_ICP_distances_Uniform[:,ICP_order]
D_WTRR_ICPU = 0.5 * (WTRR_ICP_distances_Uniform + np.transpose(WTRR_ICP_distances_Uniform))

WTRR_ICP_distances_Resize = np.loadtxt(open(result_dir /"WTRR_ICP_distances_Resize.csv", "rb"), delimiter=",")
WTRR_ICP_distances_Resize[:,:] = WTRR_ICP_distances_Resize[ICP_order,:]
WTRR_ICP_distances_Resize[:,:] = WTRR_ICP_distances_Resize[:,ICP_order]
D_WTRR_ICPS = 0.5 * (WTRR_ICP_distances_Resize + np.transpose(WTRR_ICP_distances_Resize))

WTRR_ICP_distances_Affine = np.loadtxt(open(result_dir /"WTRR_ICP_distances_Affine.csv", "rb"), delimiter=",")
WTRR_ICP_distances_Affine[:,:] = WTRR_ICP_distances_Affine[ICP_order,:]
WTRR_ICP_distances_Affine[:,:] = WTRR_ICP_distances_Affine[:,ICP_order]
D_WTRR_ICPA = 0.5 * (WTRR_ICP_distances_Affine + np.transpose(WTRR_ICP_distances_Affine))

Top_ICP_distances_Rigid = np.loadtxt(open(result_dir /"Top_ICP_distances_Rigid.csv", "rb"), delimiter=",")
Top_ICP_distances_Rigid[:,:] = Top_ICP_distances_Rigid[ICP_order,:]
Top_ICP_distances_Rigid[:,:] = Top_ICP_distances_Rigid[:,ICP_order]
D_Top_ICPR = 0.5 * (Top_ICP_distances_Rigid + np.transpose(Top_ICP_distances_Rigid))

Top_ICP_distances_Uniform = np.loadtxt(open(result_dir /"Top_ICP_distances_Uniform.csv", "rb"), delimiter=",")
Top_ICP_distances_Uniform[:,:] = Top_ICP_distances_Uniform[ICP_order,:]
Top_ICP_distances_Uniform[:,:] = Top_ICP_distances_Uniform[:,ICP_order]
D_Top_ICPU = 0.5 * (Top_ICP_distances_Uniform + np.transpose(Top_ICP_distances_Uniform))

Top_ICP_distances_Resize = np.loadtxt(open(result_dir /"Top_ICP_distances_Resize.csv", "rb"), delimiter=",")
Top_ICP_distances_Resize[:,:] = Top_ICP_distances_Resize[ICP_order,:]
Top_ICP_distances_Resize[:,:] = Top_ICP_distances_Resize[:,ICP_order]
D_Top_ICPS = 0.5 * (Top_ICP_distances_Resize + np.transpose(Top_ICP_distances_Resize))

Top_ICP_distances_Affine = np.loadtxt(open(result_dir /"Top_ICP_distances_Affine.csv", "rb"), delimiter=",")
Top_ICP_distances_Affine[:,:] = Top_ICP_distances_Affine[ICP_order,:]
Top_ICP_distances_Affine[:,:] = Top_ICP_distances_Affine[:,ICP_order]
D_Top_ICPA = 0.5 * (Top_ICP_distances_Affine + np.transpose(Top_ICP_distances_Affine))


CH_ICP_Shear = np.loadtxt(open(result_dir /"CH_ICP_Shear.csv", "rb"), delimiter=",")
CH_ICP_Shear[:,:] = CH_ICP_Shear[ICP_order,:]
CH_ICP_Shear[:,:] = CH_ICP_Shear[:,ICP_order]

Top_ICP_Shear = np.loadtxt(open(result_dir /"Top_ICP_Shear.csv", "rb"), delimiter=",")
Top_ICP_Shear[:,:] = Top_ICP_Shear[ICP_order,:]
Top_ICP_Shear[:,:] = Top_ICP_Shear[:,ICP_order]

WTRR_ICP_Shear = np.loadtxt(open(result_dir /"WTRR_ICP_Shear.csv", "rb"), delimiter=",")
WTRR_ICP_Shear[:,:] = WTRR_ICP_Shear[ICP_order,:]
WTRR_ICP_Shear[:,:] = WTRR_ICP_Shear[:,ICP_order]


WTRR_PHT_distances_Affine =  np.load(result_dir / "WTRR_PHT_distances_Affine.npy")
D_WTRR_PHTA = 0.5 * (WTRR_PHT_distances_Affine + np.transpose(WTRR_PHT_distances_Affine))

CH_PHT_distances_Affine =  np.load(result_dir / "CH_PHT_distances_Affine.npy")
D_CH_PHTA = 0.5 * (CH_PHT_distances_Affine + np.transpose(CH_PHT_distances_Affine))

WTRC_HD_distances_Affine =  np.load(result_dir / "WTRR_HD_distances_Affine.npy")
D_WTRC_HDA = 0.5 * (WTRC_HD_distances_Affine + np.transpose(WTRC_HD_distances_Affine))

WTRC_PH_distances_Affine =  np.load(result_dir / "PH Results" / "WTRC_Affine" / "distances.npy")
D_WTRC_PHA = 0.5 * (WTRC_PH_distances_Affine + np.transpose(WTRC_PH_distances_Affine))


WTRC_HD_distances_Original =  np.load(result_dir / "HD Results" / "WTRC_Original" / "distances.npy")
D_WTRC_HDO = 0.5 * (WTRC_HD_distances_Original + np.transpose(WTRC_HD_distances_Original))

CH_HD_distances_Affine =  np.load(result_dir / "HD Results" / "CH_Affine" / "distances.npy")
D_CH_HDA = 0.5 * (CH_HD_distances_Affine + np.transpose(CH_HD_distances_Affine))

CH_HD_distances_Original =  np.load(result_dir / "HD Results" / "CH_Original" / "distances.npy")
D_CH_HDO = 0.5 * (CH_HD_distances_Original + np.transpose(CH_HD_distances_Original))

CH_MHD_distances_Affine =  np.load(result_dir / "MHD Results" / "CH_Affine" / "distances.npy")
D_CH_MHDA = 0.5 * (CH_MHD_distances_Affine + np.transpose(CH_MHD_distances_Affine))

CH_MHD_distances_Original =  np.load(result_dir / "MHD Results" / "CH_Original" / "distances.npy")
D_CH_MHDO = 0.5 * (CH_MHD_distances_Original + np.transpose(CH_MHD_distances_Original))






#=========================================================================================================================
    # plot the vertices as a scatter plot
#========================================================================================================================= 
def plot(vertices):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    x =vertices[:, 0]
    y =vertices[:, 1]
    z =vertices[:, 2]



    ax.scatter(x, y, z, c='r', marker='o')

    ax.set_xlabel('X Label')
    ax.set_ylabel('Y Label')
    ax.set_zlabel('Z Label')
    
    return
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================      
    
#=========================================================================================================================
#displays a color plot of the given matrix
#=========================================================================================================================  
def plot_matrix(mat, numlabels=5, save_name='distance_Matrix_WTR.png'):
    import seaborn as sns; sns.set()
    
    
    
    labels = ['fortA', 'fortB','fortC', 'magnA', 'magnB','magnC', 'fuliA','fuliB','fuliC', 'coniA','coniB','coniC', 
              'scanA','scanB', 'septA','septB', 'paliA','paliB','paliC', 'parvA','parvB','parvC', 'psittA', 'psittB', 'psittC',
              'olivA', 'olivB', 'olivC', 'fuscA', 'fuscB', 'crassA', 'crassB', 'crassC']
    
   
    fig, ax = plt.subplots(figsize=(11.6,10))         # Sample figsize in inches
    sns.heatmap(mat, xticklabels=labels, yticklabels=labels, cmap="winter_r", ax=ax)
    
    plt.savefig(result_dir / save_name, dpi=300)
    
    plt.show()
    
    return
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================  

#=========================================================================================================================
# Applies multidimensional scaling to visualize the distance matrix
#=========================================================================================================================  
def make_MDS(distance_matrix, save_name='MDS_space_WTR.png'):
    from sklearn.manifold import MDS
    
    
    
    model = MDS(n_components=2, dissimilarity='precomputed', random_state=1)
    out = model.fit_transform(distance_matrix)    
    
    species = np.array([1,1,1, 2,2,2, 3,3,3, 4,4,4, 5,5, 6,6, 7,7,7, 8,8,8, 9,9,9, 10,10,10, 11,11, 12,12,12 ])
    
    #species = np.array([1,1,1, 1,1,1, 1,1,1, 1,1,1, 1,1, 1,1, 2,2,2, 2,2,2, 2,2,2, 3,3,3, 3,3, 4,4,4 ])
    
    fig, ax = plt.subplots(figsize=(11.6,11.6))
    
    colorize = dict(c=species, cmap=plt.cm.get_cmap('rainbow', 12)) 
    
    plt.scatter(out[:, 0], out[:, 1], **colorize, s=75)
    
    ax.tick_params(axis='both', which='major', labelsize=14)
    ax.tick_params(axis='both', which='minor', labelsize=14)
    plt.axis('equal')
    
    
    plt.savefig(result_dir / save_name, dpi=300)
    plt.show()
    
    
    return
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#========================================================================================================================= 
    

#=========================================================================================================================
# Applies multidimensional scaling to visualize the distance matrix
#=========================================================================================================================  
def make_species_legend():
   
    species = np.array([1,2, 3,4, 5,6, 7,8, 9,10,11,12 ])
    #species = np.array([1,2, 3,4])
    
    fig, ax = plt.subplots(figsize=(4,10))
    
    colorize = dict(c=species, cmap=plt.cm.get_cmap('rainbow', 12)) 
    
    x = 2*np.ones(species.shape[0]); y = (2.5*np.arange(species.shape[0])/species.shape[0]) - 1.5
    plt.axis('off')
    plt.scatter(x, y, s=100, **colorize)
    
    labels = ['  G.Fortis','  G.Magnirostris', '  G.Fuliginosa','  G.Conirostris', '  G.Scandens', '  G.Septentrionalist', 
              '  C.Palidus', '  C.Parvulus', '  C.Psittacula',  '  C2.Olivcea', '  C2.Fusca', '  P.Crassirostris']
    
    #labels = [' Geospiza', ' Camarhynchus',  ' Certhidea', ' Platyspiza']
    

    for i, txt in enumerate(labels):
        plt.annotate(txt, (x[i], y[i]), fontsize=16, fontstyle='oblique')
    
    plt.savefig(result_dir / 'species_legend.png', dpi=300)
    plt.show()
    
    
    
    return
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#========================================================================================================================= 

#=========================================================================================================================
# Applies multidimensional scaling to visualize the distance matrix
#=========================================================================================================================  
def make_MDS0(distance_matrix):
    from sklearn.manifold import MDS
    
    
      
    model = MDS(n_components=3, dissimilarity='precomputed', random_state=1)
    out = model.fit_transform(distance_matrix)    
    
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    x, y, z =out[:, 0], out[:, 1], out[:, 2]
   


    ax.scatter(x, y, z, c='r', marker='o')

    ax.set_xlabel('X Label')
    ax.set_ylabel('Y Label')
    ax.set_zlabel('Z Label')

    plt.show()
    
    
    np.savetxt(result_dir / "mds.npy", out)
    return
    
    
    species = np.array([1,1,1, 2,2,2, 3,3,3, 4,4,4, 5,5, 6,6, 7,7,7, 8,8,8, 9,9,9, 10,10,10, 11,11, 12,12,12 ])
    
    colorize = dict(c=species, cmap=plt.cm.get_cmap('rainbow', 12)) 
    
    plt.scatter(out[:, 0], out[:, 1], out[:, 2], **colorize)
   # plt.axis('equal')
    plt.show()
    
    species = np.array([1,2, 3,4, 5,6, 7,8, 9,10,11,12 ])
    
    colorize = dict(c=species, cmap=plt.cm.get_cmap('rainbow', 12)) 
    
    x = 2*np.ones(species.shape[0]); y = (2.5*np.arange(species.shape[0])/species.shape[0]) - 1.5
    plt.scatter(x, y, **colorize)
    
    labels = ['G.fortis','G.magniro', 'G.fulig','G.coniro', 'G.scand', 'G.septen', 'C.pali','C.parv', 'C.psitt',  'C2.oliv', 'C2.fusca', 
              'P.crassiro']
    

    for i, txt in enumerate(labels):
        plt.annotate(txt, (x[i], y[i]))
    
    plt.axis('equal')
    
    
    plt.show()

    return

#=========================================================================================================================
# finds the asymetrics part of the matrix, divided by the symemtric part. Zero if symmetric element is zero
#=========================================================================================================================  
def asymetric_part(mat):
    
    sym = 0.5*(mat + mat.transpose())
    
    asym = np.abs(mat - mat.transpose())
    
    rel = np.divide(asym, sym, out=np.zeros_like(asym), where=sym!=0)
    
    plot_matrix(rel)
    
    return rel.max()
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================  
#c = np.divide(a, b, out=np.zeros_like(a), where=b!=0)

#names of beaks
beak1 =  "G.FortisA"
beak2 =  "G.FortisB"
beak3 =  "G.FortisC"
beak4 = "G.MagnirostrisA"
beak5 =  "G.MagnirostrisB"
beak6 =  "G.MagnirostrisC"
beak7 =  "G.FuliginosaA"
beak8 = "G.FuliginosaB"
beak9 = "G.FuliginosaC"
beak10 =  "G.ConirostrisA"
beak11 =  "G.ConirostrisB"
beak12 =  "G.ConirostrisC"
beak13 = "G.ScandensA"
beak14 =   "G.ScandensB"
beak15 =  "G.SeptentrionalistA"
beak16 =  "G.SeptentrionalistB"
beak17 =  "C.PallidusA"
beak18 =  "C.PallidusB"
beak19 = "C.PallidusC"
beak20 =  "C.ParvulusA"
beak21 =  "C.ParvulusB"
beak22 =  "C.ParvulusC"
beak23 =  "C.PsittaculaA"
beak24 = "C.PsittaculaB"
beak25 =  "C.PsittaculaC"
beak26 = "C2.OlivaceaA"
beak27 =  "C2.OlivaceaB"
beak28 =  "C2.OlivaceaC"
beak29 =  "C2.FuscaA"
beak30 =  "C2.FuscaB"
beak31 =  "P.CrassirostrisA"
beak32 =  "P.CrassirostrisB"
beak33 =  "P.CrassirostrisC"
names = np.array([beak1, beak2, beak3, beak4, beak5, beak6, beak7, beak8, beak9, 
                  beak10, beak11, beak12, beak13, beak14, beak15, beak16, beak17, beak18,
                  beak19, beak20, beak21, beak22, beak23, beak24, beak25, beak26, beak27,
                  beak28, beak29, beak30, beak31, beak32, beak33]) 



def path_leaf(path):
    import ntpath
    head, tail = ntpath.split(path)
    return tail or ntpath.basename(head)


def get_transformations(root_dir, name):
# Return a list of regular files only, not directories
    file_list = [f for f in root_dir.glob('**/*') if f.is_file()]
    
    transformations = np.zeros((33,33, 4, 4))
    
    for f in file_list: 
        for i, beak1 in enumerate(names):
            for j, beak2 in enumerate(names):
                if i==j:
                    transformations[i, j] = np.eye(4)
                
                elif path_leaf(f).startswith(beak1) and (beak2 in path_leaf(f)):
                    transformations[i, j] = np.loadtxt(f)
                    
    np.save(result_dir/ ('transformations_'+ name + '.npy'), transformations)
    return transformations


def collapse_meshes(root_dir, transformations):
    '''
    root_dir contrains the meshes to be transformed. 
    transformations is a vector of transformation matrices
    '''
    #need to find the right mesh to transform based on the order of the name
    
    file_list = [f for f in root_dir.glob('**/*') if f.is_file()]
    
    num_beaks = len(names)
    
    for i in range(num_beaks):
        for j in range(num_beaks):
            beak1 = file_list[i]
            beak2 = names[j]
            
            if (beak2 in path_leaf(beak1)):
                transformation = transformations[j]
                matrix = np.transpose(transformation[:3, :3])
                translation =  transformation[:3, 3]
                AM.apply_transformation(beak1, matrix, translation)
                
    return
    
    
    
    
    
beak1 =  mesh_dir / "G.FortisA_WTR.stl"
beak2 =  mesh_dir / "G.FortisB_WTR.stl"
beak3 =  mesh_dir / "G.FortisC_WTR.stl"
beak4 = mesh_dir / "G.MagnirostrisA_WTR.stl"
beak5 = mesh_dir / "G.MagnirostrisB_WTR.stl"
beak6 = mesh_dir / "G.MagnirostrisC_WTR.stl"
beak7 = mesh_dir /"G.FuliginosaA_WTR.stl"
beak8 = mesh_dir /"G.FuliginosaB_WTR.stl"
beak9 = mesh_dir /"G.FuliginosaC_WTR.stl"
beak10 =  mesh_dir / "G.ScandensA_WTR.stl"
beak11 =  mesh_dir / "G.ScandensB_WTR.stl"
beak12 =  mesh_dir / "G.SeptentrionalistA_WTR.stl"
beak13 =  mesh_dir / "G.SeptentrionalistB_WTR.stl"
beak14 = mesh_dir / "C.PallidusA_WTR.stl"
beak15 = mesh_dir / "C.PallidusB_WTR.stl"
beak16 = mesh_dir / "C.PallidusC_WTR.stl"
beak17 = mesh_dir / "C.PsittaculaA_WTR.stl"
beak18 = mesh_dir / "C.PsittaculaB_WTR.stl"
beak19 = mesh_dir / "C.PsittaculaC_WTR.stl"
beak20 = mesh_dir / "C.ParvulusA_WTR.stl"
beak21 = mesh_dir / "C.ParvulusB_WTR.stl"
beak22 = mesh_dir / "C.ParvulusC_WTR.stl"
beak23 = mesh_dir / "C2.FuscaA_WTR.stl"
beak24 = mesh_dir / "C2.FuscaB_WTR.stl"
beak25 = mesh_dir / "C2.OlivaceaA_WTR.stl"
beak26 = mesh_dir / "C2.OlivaceaB_WTR.stl"
beak27 = mesh_dir / "C2.OlivaceaC_WTR.stl"
beak28 =  mesh_dir / "G.ConirostrisA_WTR.stl"
beak29 =  mesh_dir / "G.ConirostrisB_WTR.stl"
beak30 =  mesh_dir / "G.ConirostrisC_WTR.stl"
beak31 =  mesh_dir / "P.CrassirostrisA_WTR.stl"
beak32 =  mesh_dir / "P.CrassirostrisB_WTR.stl"
beak33 =  mesh_dir / "P.CrassirostrisC_WTR.stl"
beaks_wtr = np.array([beak1, beak2, beak3, beak4, beak5, beak6, beak7, beak8, beak9, 
                  beak10, beak11, beak12, beak13, beak14, beak15, beak16, beak17, beak18,
                  beak19, beak20, beak21, beak22, beak23, beak24, beak25, beak26, beak27,
                  beak28, beak29, beak30, beak31, beak32, beak33]) 
        
        
beak1 =  mesh_dir / "G.FortisA.stl"
beak2 =  mesh_dir / "G.FortisB.stl"
beak3 =  mesh_dir / "G.FortisC.stl"
beak4 = mesh_dir / "G.MagnirostrisA.stl"
beak5 = mesh_dir / "G.MagnirostrisB.stl"
beak6 = mesh_dir / "G.MagnirostrisC.stl"
beak7 = mesh_dir /"G.FuliginosaA.stl"
beak8 = mesh_dir /"G.FuliginosaB.stl"
beak9 = mesh_dir /"G.FuliginosaC.stl"
beak10 =  mesh_dir / "G.ScandensA.stl"
beak11 =  mesh_dir / "G.ScandensB.stl"
beak12 =  mesh_dir / "G.SeptentrionalistA.stl"
beak13 =  mesh_dir / "G.SeptentrionalistB.stl"
beak14 = mesh_dir / "C.PallidusA.stl"
beak15 = mesh_dir / "C.PallidusB.stl"
beak16 = mesh_dir / "C.PallidusC.stl"
beak17 = mesh_dir / "C.PsittaculaA.stl"
beak18 = mesh_dir / "C.PsittaculaB.stl"
beak19 = mesh_dir / "C.PsittaculaC.stl"
beak20 = mesh_dir / "C.ParvulusA.stl"
beak21 = mesh_dir / "C.ParvulusB.stl"
beak22 = mesh_dir / "C.ParvulusC.stl"
beak23 = mesh_dir / "C2.FuscaA.stl"
beak24 = mesh_dir / "C2.FuscaB.stl"
beak25 = mesh_dir / "C2.OlivaceaA.stl"
beak26 = mesh_dir / "C2.OlivaceaB.stl"
beak27 = mesh_dir / "C2.OlivaceaC.stl"
beak28 =  mesh_dir / "G.ConirostrisA.stl"
beak29 =  mesh_dir / "G.ConirostrisB.stl"
beak30 =  mesh_dir / "G.ConirostrisC.stl"
beak31 =  mesh_dir / "P.CrassirostrisA.stl"
beak32 =  mesh_dir / "P.CrassirostrisB.stl"
beak33 =  mesh_dir / "P.CrassirostrisC.stl"
beaks_t = np.array([beak1, beak2, beak3, beak4, beak5, beak6, beak7, beak8, beak9, 
                  beak10, beak11, beak12, beak13, beak14, beak15, beak16, beak17, beak18,
                  beak19, beak20, beak21, beak22, beak23, beak24, beak25, beak26, beak27,
                  beak28, beak29, beak30, beak31, beak32, beak33]) 




matrix1 =  matrix_dir / "C.PallidusA_G.FortisA.txt"
matrix2 =  matrix_dir / "C.PallidusA_G.FortisB.txt"
matrix3 =  matrix_dir / "C.PallidusA_G.FortisC.txt"
matrix4 =  matrix_dir / "C.PallidusA_G.MagnirostrisA.txt"
matrix5 =  matrix_dir / "C.PallidusA_G.MagnirostrisB.txt"
matrix6 =  matrix_dir / "C.PallidusA_G.MagnirostrisC.txt"
matrix7 =  matrix_dir / "C.PallidusA_G.FuliginosaA.txt"
matrix8 =  matrix_dir / "C.PallidusA_G.FuliginosaB.txt"
matrix9 =  matrix_dir /  "C.PallidusA_G.FuliginosaC.txt"
matrix10 =  matrix_dir / "C.PallidusA_G.ScandensA.txt"
matrix11 =  matrix_dir / "C.PallidusA_G.ScandensB.txt"
matrix12 =  matrix_dir / "C.PallidusA_G.SeptentrionalistA.txt"
matrix13 =  matrix_dir / "C.PallidusA_G.SeptentrionalistB.txt"
matrix14 =  matrix_dir / "C.PallidusA_C.PallidusA.txt"
matrix15 =  matrix_dir / "C.PallidusA_C.PallidusB.txt"
matrix16 =  matrix_dir / "C.PallidusA_C.PallidusC.txt"
matrix17 =  matrix_dir / "C.PallidusA_C.PsittaculaA.txt"
matrix18 =  matrix_dir /  "C.PallidusA_C.PsittaculaB.txt"
matrix19 =  matrix_dir / "C.PallidusA_C.PsittaculaC.txt"
matrix20 =  matrix_dir / "C.PallidusA_C.ParvulusA.txt"
matrix21 =  matrix_dir / "C.PallidusA_C.ParvulusB.txt"
matrix22 =  matrix_dir / "C.PallidusA_C.ParvulusC.txt"
matrix23 =  matrix_dir / "C.PallidusA_C2.FuscaA.txt"
matrix24 =  matrix_dir / "C.PallidusA_C2.FuscaB.txt"
matrix25 =  matrix_dir / "C.PallidusA_C2.OlivaceaA.txt"
matrix26 =  matrix_dir / "C.PallidusA_C2.OlivaceaB.txt"
matrix27 =  matrix_dir /  "C.PallidusA_C2.OlivaceaC.txt"
matrix28 =  matrix_dir / "C.PallidusA_G.ConirostrisA.txt"
matrix29 =  matrix_dir / "C.PallidusA_G.ConirostrisB.txt"
matrix30 =  matrix_dir / "C.PallidusA_G.ConirostrisC.txt"
matrix31 =  matrix_dir / "C.PallidusA_P.CrassirostrisA.txt"
matrix32 =  matrix_dir / "C.PallidusA_P.CrassirostrisB.txt"
matrix33 =  matrix_dir / "C.PallidusA_P.CrassirostrisC.txt"
matrices = np.array([matrix1, matrix2, matrix3, matrix4, matrix5, matrix6, matrix7, matrix8, matrix9, 
                  matrix10, matrix11, matrix12, matrix13, matrix14, matrix15, matrix16, matrix17, matrix18,
                  matrix19, matrix20, matrix21, matrix22, matrix23, matrix24, matrix25, matrix26, matrix27,
                  matrix28, matrix29, matrix30, matrix31, matrix32, matrix33]) 

'''
for i,  beak in enumerate(beaks_t):
    transformation = np.loadtxt(matrices[i])
    matrix = np.transpose(la.inv(transformation[:3, :3]))
    translation = - transformation[:3, 3]
    AM.apply_transformation(beak, matrix, translation)
'''


'''
result_dir1 = result_dir / "Run with CH Small Cut Not Aligned"
result_dir2 = result_dir / "Run with WTR Aligned"
result_dir3 = result_dir / "Run with CH Small Cut Aligned"

distances1 = np.load(result_dir1 / "distances.npy")

transformations1 = np.load(result_dir1 / "transformations.npy")


distances2 = np.load(result_dir2 / "distances.npy")
D2 = 0.5 * (distances2 + np.transpose(distances2))

transformations2 = np.load(result_dir2 / "transformations_scale_opt.npy")

distances3 = np.load(result_dir3 / "distances_affine_opt.npy")

transformations3 = np.load(result_dir3 / "transformations_affine_opt.npy")

shears1 = np.zeros_like(distances1) 
for i in range(distances1.shape[0]):
    for j in range(distances1.shape[0]):
        sl3 = transformations1[i,j]
        shears1[i,j] = max(np.abs([sl3[0,1],sl3[0,2],sl3[1,2]]))

shearsA3 = np.zeros_like(distances3) 
shearsB3 = np.zeros_like(distances3) 
for i in range(distances3.shape[0]):
    for j in range(distances1.shape[0]):
        sl3 = transformations3[i,j]
        
        shearsA3[i,j] = max(np.abs([sl3[0,1],sl3[0,2],sl3[1,2]]))
        shearsB3[i,j] = max(np.abs([sl3[0,1],sl3[0,2],sl3[1,2]]))/min(np.abs([sl3[0,0],sl3[1,1],sl3[2,2]]))
        
dets3 = np.zeros_like(distances3) 
for i in range(distances3.shape[0]):
    for j in range(distances1.shape[0]):
        sl3 = transformations3[i,j]
        dets3[i,j] = np.abs(la.det(sl3))
'''















